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We analyze stability of superfluid currents in a system of strongly interacting ultra-cold atoms in 
an optical lattice. We show that such a system undergoes a dynamic, irreversible phase transition at 
a critical phase gradient that depends on the interaction strength between atoms. At commensurate 
filling, the phase boundary continuously interpolates between the classical modulation instability of 
a weakly interacting condensate and the equilibrium quantum phase transition into a Mott insulator 
state at which the critical current vanishes. We argue that quantum fluctuations smear the transition 
boundary in low dimensional systems. Finally we discuss the implications to realistic experiments. 



The quantum phase transition from superfluid 
(SF) to Mott insulator (IN) ^ is an important 
paradigm of strongly correlation physics. Recently 
this transition was demonstrated in spectacular ex- 
periments involving ultra cold atoms in an optical 
lattice 2, 3J. An important feature of these sys- 
tems is that they can be essentially isolated from 
the environment, opening unique possibilities to 
study, not only the equilibrium phase diagram, but 
also quantum dynamics very far from equilibrium 
li, 1111, [3,11, la, llll- In particular, it is now possi- 
ble to explore a new class of phenomena involving 
non-equilibrium dynamics in the vicinity of quan- 
tum phase transitions. 

In this letter we analyze the stability of super- 
fluid current flow in a system of strongly interact- 
ing bosons in a lattice. We show that at the mean 
field level such a system, undergoes a dynamic 
phase transition, associated with irreversible de- 
cay of the superfluid flow at a critical momentum 
that depends on the interaction strength. We ar- 
gue that quantum phase fluctuations play an im- 
portant role near the phase boundary. In systems 
of lower dimensionality, they broaden the transi- 
tion significantly as was indeed observed in recent 
experiments Tuj and numerical simulations lllLll2l | 
in one dimensional systems. In three dimensions, 
by contrast, we predict that the current decay rate 
exhibits a sharp discontinuity at the mean field 
transition. We propose to test this prediction with 
transport experiments in three dimensional optical 
lattices . 

A well studied effect, closely related to our dis- 
cussion, is the modulational instability of weakly 
interacting bosons on a lattice |g. It was exper- 
imentally demonstrated d ^3 ^^^^ ^^ch a con- 
densate undergoes a dynamical localization transi- 
tion, involving onset of chaos, when the phase gra- 
dient (or condensate momentum) associated with 



the flow exceeds 7r/2 per lattice unit. The dynamic 
phase transition described in this letter interpo- 
lates continuously between the classical instability 
at condensate momentum 7r/2 and the quantum 
phase transition into the Mott state at zero cur- 
rent, thereby establishing a natural connection be- 
tween the two transitions. 



Dependence of the critical momentum on the in- 
teraction can be understood as follows. The su- 
perfluid current / associated with a condensate 
moving within the lowest Bloch band is I{p) ^ 
Ps sin(p), where p is the (quasi)momentum of the 
condensate measured in the units of inverse lat- 
tice constant and ps is the superfluid density. In 
a weakly interacting condensate ps is independent 
of p. Thus, the current increases with p up to a 
maximal value at = 7r/2. Beyond this point, 
the effective mass changes sign and any further in- 
crease in p results in decrease of the current, ren- 
dering the superfluid unstable At strong in- 
teractions p itself is a function of the effective mass 
and thus also of p. Specifically p decreases as p is 
increased, such that the maximum of I{p) occurs 
in general a,t Pc < Tr/2. In particular, close to the 
Mott insulator, a slight increase in effective mass, 
and thus also in p leads to vanishing of p. There- 
fore Pc tends to zero toward the Mott transition. 
These considerations include the effect of quantum 
depletion at the mean field level, which is effective 
in all dimensions, and lead to a stability phase dia- 
gram (for commensurate filling) depicted in Fig. ^ 
Smearing of the transition lines by quantum phase 
slips beyond mean field theory, is effective only in 
lower dimensionality. 

Static and dynamic properties of condensates in 
optical lattices are described by the Bose Hubbard 
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FIG. 1: Stability phase diagram in the plane of phase 
twist per bond versus dimensionless interaction for fill- 
ing of A'^ = 1 particle per site found from numerical 
solution of time dependent Gutzwiller equations JSJl. 



model (BHM) 



where J is the hopping amplitude between the 
nearest neighbors (ij), U is the on-site repulsive 
interaction, and = a|ai is the number operator. 
We denote the average number of bosons per site 
by N. Provided the system is deep in the super- 
fluid phase {JN 3> U), the dynamics can be well 
approximated by the discrete GP equations 
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where '0i = (ij) is the matter field and the set 
O contains the nearest neighbors of site i. Linear 
mode analysis around the stationary current carry- 
ing solutions ijji = \/N exoU/pXi), yields the onset 
of instability at p = Tr/2 |3, l3- 

Close to the Mott transition, increased quantum 
fluctuations invalidate the GP description. How- 
ever, one can still use semiclassical order parame- 
ter dynamics if one coarse grains the system into 
blocks of roughly a coherence length ^, which is re- 
lated to the superfluid density by a standard scal- 
ing form 16]. As in the weak coupling theory, the 
current should become unstable when the phase 
change per unit cell exceeds 7r/2, with the unit cell 



now of order ^. Since the coherence length, ^, di- 
verges at the transition, we expect Pc to vanish as 
1/^, as we indeed find below. The diverging length 
scale facilitates a continuum description of the dy- 
namics close to the Mott transition, in the form of a 
time dependent Ginzburg-Landau equation Q| 
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Within the mean-field approximation, and in the 
limit of large average occupation 3> 1, C — 
2d{l — u), where u = U / {8JNd) is the dimension- 
less interaction constant. If N is not too large then 
Eq.© still holds but the expressions for ^ and u 
are more complicated Q. 

We choose the zero of energy at the state with 
integer filling. Then the deviation from commen- 
surate density is given in terms of the superfluid 
order parameter ip hy Sn = Cd{ip*ip — V''0*)/2«, 
with Cd — u^^(2(i)~^/^ and d being the dimen- 
sionality of the system. Note that / d'^xSn is a 
constant of motion under (O • The supercurrent is 
similarly given by / = Cd{'>p*'^ip — i)V ip*) / 2i . 
The equation of motion ^ admits uniform so- 
pe^P'^+^t'^ with p 
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p^ ~ p^ 



which are characterized by a phase gradient p and 
a relative density 



6n = CdKC 



p ~p 
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In particular, at commensurate filling p — Q and 
ip is time-independent. To analyze whether these 
solutions are stable we find the spectrum of small 
fluctuations around them. There are two eigen- 
modes in the superfluid regime (^~^ > 0): a sta- 
ble gapped mode, and a phase (Bogoliubov) mode 
with linear dispersion at long wave lengths. The 
spectrum of the latter, for wavevectors parallel to 
the current, reads: 



uj{k) 



2pp 



2^2 _|_ p2 2/i2 + p2 
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Here the first term is analogous to the usual 
Doppler shift and the second describes propagation 
of the sound waves in the moving reference frame. 
The onset of imaginary frequencies marking the in- 
stability occurs at p2_^(p^)2 _ 1/3^2, Combining 

this with Eq. ^ we find that for iV > 1 
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As argued before on general grounds, in the case of 
commensurate filling Sn = 0, the critical phase gra- 
dient vanishes toward the equilibrium Mott transi- 
tion (m = 1) as Pc oc 1/^ oc -s/l — u. At incommen- 
surate density there is no equilibrium Mott tran- 
sition. As a result, we do not expect the insta- 
bility to reach p ^ 0. Indeed, Pc has a minimum 
at u < 1 < oo), and diverges as u 1. The 
divergence, simply signals the breakdown of the 
continuum theory and is cutoff by the lattice. 

To interpolate between the regimes of weak and 
strong interactions we employ the Gutzwiller ap- 
proximation [isj . In this approach, the wavefunc- 
tion is assumed to be factorizable: 
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Here j is a site index and n is the site occupation. 
The ansatz ((Tjl supplemented by self-consistency 
conditions leads to equations of motion for the vari- 
ational parameters: 



where 
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For actual calculations we truncated {Tj) at five and 
ten states per site without noticeable differences in 
the results. 

Equations 0, admit uniform current carrying 
solutions. We numerically check their stability to 
slight perturbations in the equations of motion. 
We show the stability boundaries at commensurate 
filling in Fig. ^ It is evident that the dynamical 
instability at 7r/2 in the GP regime is continuously 
connected to the equilibrium (zero current) Mott 
transition. 

We perform a similar analysis at incommensu- 
rate filling (Fig. ^ . In agrement with the contin- 
uum expression © we find that the critical mo- 
mentum Pc reaches a minimum at some u < 1. At 
stronger interactions, Pc increases, and saturates 
at tt/2. 

The mean field transition discussed above ig- 
nores the possibility of current decay below the 
critical momentum due to quantum tunneling out 
of the metastable state. Such processes are expo- 
nentially suppressed by a tunneling action through 
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FIG. 2: Stability phase diagrams for different filling 
factors in a two-dimensional lattice from numerical so- 
lution of the time dependent Gutzwiller equations l|HJ. 
Away from commensurate filling the critical phase gra- 
dient reaches a minimum and climbs back to 7r/2. 



a barrier. But since the barrier vanishes at the 
mean field instability, they can potentially broaden 
the transition rendering the phase diagrams of 
Figs, n] and |21 meaningless. This problem is ad- 
dressed in full detail in Ref. within the general 
framework of Here we one important result 

for the decay rate close to the equilibrium SF-IN 
transition, i.e. at small critical current. To obtain 
the tunnelling action close to the critical current we 
expand the GL action associated with Q around 
the metastable solution to cubic order in the fluc- 
tuations. Then we use a scaling approach similar 
to the one introduced in Ref. [l8| in the context 
of spinodal decomposition. Since the instability 
in (0) arises at k — s- the barrier vanishes only 
for a tunneling instanton of diverging (space-time) 
volume. By dimensional analysis we find Vinst oc 
{pc — pY^^I"^ ^ while the energy density of the bar- 
rier E\, (X (pc — vf" ■ Consequently the instanton 
action is Sinst ^ Et x V^nst oc Bd{pc - p)^'^"''. 
Note that this scaling is different from that derived 
in Ref. because the field if) is complex. Thus 
in one and two dimensions the tunneling action 
vanishes continuously toward the mean field tran- 
sition. Then we expect the transition to become a 
wide crossover {Bd is calculated in Q| and found 
to be of order unity) . This agrees with recent ex- 
periments and numerical simulations TTl in 
one dimensional traps. 
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In three dimensions, by contrast, the action of 
"critical" instantons diverges because of their di- 
verging volume. Current would then rather decay 
via non critical instantons, i.e. ones of finite size 
that feel a finite energy barrier, and thus cost finite 
action. We therefore predict that in three dimen- 
sions, the decay rate will exhibit a discontinuity 
at the mean field transition. A variational calcula- 
tion 01 yields a rate T oc e~'^ '^ in the vicinity of 
the transition. In this sense the mean-field phase 
diagram (Fig. ^) is well defined in three dimen- 
sions. 

In realistic experimental situations condensates 
are confined in harmonic traps which leads to a 
non uniform density distribution in the form of do- 
mains with different N. In the weakly interacting 
regime the critical momentum is Pc — 7r/2, i.e. it is 
insensitive to the spatial density variation induced 
by the harmonic confinement. By contrast, in the 
regime of strong interactions the position of the dy- 
namical instability strongly depends on the filling 
factor N that is directly affected by the density dis- 
tribution. In particular, the motion first becomes 
unstable for the smallest integer filling iV = 1. 

In Fig. O we plot the time evolution of the 
condensate momentum (computed within the 
Gutzwiller approximation) in a trap, for two dif- 
ferent filling factors. The center of mass motion 
becomes unstable at approximately the same in- 
teraction strength in both cases. But while at 
smaller filling the condensate motion rapidly be- 
comes chaotic as in the uniform case, damping of 
oscillations at larger filling occurs much more grad- 
ually. These results can be understood by noting 
that if the phase gradient in the condensate ex- 
ceeds the critical value corresponding to TV = 1 
these domains become unstable triggering the de- 
cay of current. However, when there is high filling 
of the central sites the overall weight of domains 
with = 1, is small. Thus the effect of the insta- 
bility on the total current is reduced. 

An important experimental manifestation of 
these results is the inherently irreversible nature 
of the phase transition at finite currents. Consider 
a situation in which a moving condensate is first 
prepared on a weak lattice. Then, the depth of the 
periodic potential is increased adiabatically [l^ . 
which corresponds to moving along a horizontal 
line in the parameter space of Fig. Finally, 
the lattice depth is slowly decreased back to its 
original state and the visibility of the interference 
fringes compared to their initial value. If in this 
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FIG. 3: Time dependence of the condensate momen- 
tum in a two-dimensional harmonic trap with different 
filling factors per central site. The simulated system is 
a lattice of dimensions 120 x 60 with global trapping 
potential V{jx,jy) = 0.01(j^ + Jy)- We set the hop- 
ping amplitude J = 1 while increasing the interaction 
linearly in time: U{t) = O.Oli. 



sequence we pass the instability, then the current 
will decay into incoherent excitations and heat the 
condensate. This will result in total loss of current 
and reduced visibility of the interference fringes at 
the end of the cycle. Such experiments could be 
used to probe the nonequilibrium phase diagram 
(Fig. ^ and to determine the position of the equi- 
librium Mott transition. 
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